{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RRTM Tropical (TRP) standard atmosphere"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This notebook reads in the TRP standard atmosphere, then allows the user to scale the profiles for each of the seven gases (H2O, CO2, O3, N2O, CO, CH4, and O2).  It is also possible to scale the profiles by a factor of zero to eliminate that gas from the RRTM calculation.  It also allows the user to specify aerosols or cloud properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 1) DECISIONS REGARDING MODEL ATMOSPHERE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ......................................................................\n",
    "# Create a unique file extension for your input and output files (e.g., '2xCO2')\n",
    "fileExt     = 'trp_sw' \n",
    "\n",
    "# ......................................................................\n",
    "# Shortwave (sw) or longwave (lw) calculation?\n",
    "type = 'sw'        # Choose either 'sw' or 'lw'.\n",
    "aer  = 'n'         # Choose either 'n' or 'y'.\n",
    "cld  = 'n'         # Currently one must choose either aerosols OR clouds, but not both.\n",
    "\n",
    "# ......................................................................\n",
    "# Solar zenith angle\n",
    "sza  = 30.0\n",
    "\n",
    "# ......................................................................\n",
    "# Scalings for the atmospheric gases?   (Set all to 1.0 for NO SCALING.)\n",
    "h2o_scale = 1.0\n",
    "co2_scale = 1.0\n",
    "o3_scale  = 1.0\n",
    "n2o_scale = 1.0\n",
    "co_scale  = 1.0\n",
    "ch4_scale = 1.0\n",
    "o2_scale  = 1.0\n",
    "\n",
    "# ......................................................................\n",
    "# Include aerosols?   (If 'y', then you must specify information for the IN_AER_RRTM file below.)\n",
    "if aer=='y':\n",
    "    aerFilename = 'IN_AER_RRTM'\n",
    "    fa = open(aerFilename,'w')\n",
    "    fa.write('    1\\n')\n",
    "    fa.write('    4    0    0    0   2.184    1.00    0.00\\n')  # Last 3 values specify the wavelength dependence; see next cell.\n",
    "    fa.write('    1 0.0013\\n')                                  # Layer 1, Aerosol optical depth at 1 micron\n",
    "    fa.write('    2 0.0037\\n')                                  # Layer 2, Aerosol optical depth at 1 micron\n",
    "    fa.write('    3 0.0037\\n')                                  # Layer 3, Aerosol optical depth at 1 micron\n",
    "    fa.write('    4 0.0037\\n')                                  # Layer 4, Aerosol optical depth at 1 micron\n",
    "    fa.write('0.75\\n')                                          # Single scattering albedo\n",
    "    fa.write('0.700\\n')                                         # Asymmetry factor (g)\n",
    "    fa.close()\n",
    "\n",
    "# ......................................................................\n",
    "# Include clouds?     (If 'y', then you must specify information for the IN_CLD_RRTM file below.)\n",
    "if cld=='y':\n",
    "    cldFilename = 'IN_CLD_RRTM'\n",
    "    fc = open(cldFilename,'w')\n",
    "    fc.write('    2    3    1\\n')                                           # Standard specification for both sw and lw.\n",
    "    fc.write('   11       1.0     6.515       0.5      65.0      10.0\\n')   # LAY, CLDFRAC, CWP, FRACICE, EFFSIZEICE, EFFSIZELIQ\n",
    "    fc.write('%\\n')\n",
    "    fc.close()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 2) READ THE MODEL ATMOSPHERE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read model atmosphere, two lines at a time.\n",
    "data = loadtxt('trp.dat')\n",
    "nalt = len(data)\n",
    "\n",
    "# Extract the profile data.\n",
    "Z    = data[:,1]\n",
    "P    = data[:,2]\n",
    "T    = data[:,3]\n",
    "\n",
    "h2o  = data[:,6]\n",
    "co2  = data[:,7]\n",
    "o3   = data[:,8]\n",
    "n2o  = data[:,9]\n",
    "co   = data[:,10]\n",
    "ch4  = data[:,11]\n",
    "o2   = data[:,12]\n",
    "\n",
    "# Units for the different profiles.\n",
    "un1  = 'AA'\n",
    "un2  = 'AAAAAAA'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 3) APPLY THE SCALE FACTORS FOR EACH OF THE SEVEN GASES."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the scale factors.\n",
    "h2o       = h2o * h2o_scale\n",
    "co2       = co2 * co2_scale\n",
    "o3        = o3  * o3_scale\n",
    "n2o       = n2o * n2o_scale\n",
    "co        = co  * co_scale\n",
    "ch4       = ch4 * ch4_scale\n",
    "o2        = o2  * o2_scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 4) CREATE A NEW MODEL ATMOSPHERE (INPUT_RRTM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create A NEW text file.\n",
    "\n",
    "# NOTE: Any existing input file is deleted !!\n",
    "!rm INPUT_RRTM\n",
    "\n",
    "# Create new input file.\n",
    "newFilename = 'INPUT_RRTM'\n",
    "fn = open(newFilename,'w')\n",
    "\n",
    "if type=='sw':\n",
    "    fn.write('INPUT_RRTM_SW FOR TRP ATMOSPHERE\\n')\n",
    "    fn.write('0        1         2         3         4         5         6         7         8         9\\n')\n",
    "    fn.write('123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789-\\n')\n",
    "    fn.write('$ TROPICAL ATMOSPHERE\\n')\n",
    "    if aer=='y':\n",
    "        fn.write('                   0                             1                                0 1   00    0   10\\n')\n",
    "    elif cld=='y':\n",
    "        fn.write('                   0                             1                                0 1   00    1   10\\n')        \n",
    "    else:\n",
    "        fn.write('                   0                             1                                0 1   00    0   10\\n')\n",
    "    fn.write('              0%10.3f\\n' % sza)\n",
    "    fn.write('           2  00.8810.7940.7380.7270.6570.6600.6260.6550.7910.8830.9570.9580.9580.970\\n')\n",
    "    fn.write('    0        38         0    7    0    0     0.000     0.000    45.000\\n')\n",
    "    fn.write('     0.000    60.000\\n')\n",
    "    fn.write('     0.000     1.000     2.000     3.000     4.000     5.000     6.000     7.000\\n')\n",
    "    fn.write('     8.000     9.000    10.000    11.000    12.000    13.000    14.000    15.000\\n')\n",
    "    fn.write('    16.000    17.000    18.000    19.000    20.000    21.000    22.000    23.000\\n')\n",
    "    fn.write('    24.000    25.000    27.500    30.000    32.500    35.000    37.500    40.000\\n')\n",
    "    fn.write('    42.500    45.000    47.500    50.000    55.000    60.000\\n')\n",
    "    fn.write('   38 Tropical Atmosphere\\n')\n",
    "else:\n",
    "    fn.write('INPUT_RRTM_LW FOR TRP ATMOSPHERE\\n')\n",
    "    fn.write('0        1         2         3         4         5         6         7         8         9\\n')\n",
    "    fn.write('123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789-\\n')\n",
    "    fn.write('$ TROPICAL ATMOSPHERE\\n')\n",
    "    if cld=='n':\n",
    "        fn.write('                                                 1                   0            0 0    0    0\\n')\n",
    "    else:\n",
    "        fn.write('                                                 1                   0            2 0    0    1\\n')\n",
    "    fn.write(' 294.2\\n')\n",
    "    fn.write('    0        38         0    7    0    0     0.000                         0.000\\n')\n",
    "    fn.write('     0.000    60.000\\n')\n",
    "    fn.write('     0.000     1.000     2.000     3.000     4.000     5.000     6.000     7.000\\n')\n",
    "    fn.write('     8.000     9.000    10.000    11.000    12.000    13.000    14.000    15.000\\n')\n",
    "    fn.write('    16.000    17.000    18.000    19.000    20.000    21.000    22.000    23.000\\n')\n",
    "    fn.write('    24.000    25.000    27.500    30.000    32.500    35.000    37.500    40.000\\n')\n",
    "    fn.write('    42.500    45.000    47.500    50.000    55.000    60.000\\n')\n",
    "    fn.write('   38 Tropical Atmosphere\\n')\n",
    "\n",
    "# Write the NEW model atmosphere.\n",
    "# ....Remaining lines in pairs\n",
    "for alt in range(nalt):\n",
    "    fn.write('%10.3f%10.3f%10.3f     %2s   %7s\\n' % (Z[alt], P[alt], T[alt], un1, un2))\n",
    "    fn.write('%10.3E%10.3E%10.3E%10.3E%10.3E%10.3E%10.3E\\n' % (h2o[alt], co2[alt], o3[alt], n2o[alt], co[alt], ch4[alt], o2[alt]))\n",
    "\n",
    "# End of input file.\n",
    "fn.write('%%%%%')\n",
    "\n",
    "fn.close()\n",
    "\n",
    "# Save the input file using the desired extension.\n",
    "fileExt = '.' + fileExt\n",
    "!cp INPUT_RRTM INPUT_RRTM$fileExt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 5) RUN RRTM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Link to the executable file, LBLRTM\n",
    "!ln -s /home/vonw/radtran/rrtm/rrtm_sw/rrtm_sw_linux_pgi_v2.5 rrtm_sw\n",
    "!ln -s /home/vonw/radtran/rrtm/rrtm_lw/rrtm_v3.3_linux_pgi rrtm_lw\n",
    "\n",
    "# Link to essential input files, INPUT_RRTM (model atmosphere)\n",
    "if type=='sw':\n",
    "    !./rrtm_sw\n",
    "    !tail -n+6 OUTPUT_RRTM | head -n-17 > OUTPUT_RRTM$fileExt\n",
    "else:\n",
    "    !./rrtm_lw\n",
    "    !tail -n+4 OUTPUT_RRTM | head -n-21 > OUTPUT_RRTM$fileExt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 6) QUICK PLOT OF THE RESULTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the output file and graph the data.\n",
    "out   = loadtxt('OUTPUT_RRTM' + fileExt)\n",
    "\n",
    "if type=='sw':\n",
    "    p     = out[:,1]\n",
    "    Fu    = out[:,2]\n",
    "    Fd_df = out[:,3]\n",
    "    Fd_dr = out[:,4]\n",
    "    Fd    = out[:,5]\n",
    "    Fnet  = out[:,6]\n",
    "    hr    = out[:,7]\n",
    "    \n",
    "    plot(Fu,p,Fd_df,p,Fd_dr,p,Fd,p,Fnet,p)\n",
    "    axis([0, 1400, 0, 1013])\n",
    "    ax = gca()\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    xlabel('Flux (W m-2)')\n",
    "    ylabel('Pressure (mb)')\n",
    "    title('Shortwave Fluxes')\n",
    "    legend(('F_up','F_down_diffuse','F_down_direct','F_down','F_net'),loc='best')\n",
    "else:\n",
    "    p     = out[:,1]\n",
    "    Fu    = out[:,2]\n",
    "    Fd    = out[:,3]\n",
    "    Fnet  = out[:,4]\n",
    "    hr    = out[:,5]\n",
    "    \n",
    "    plot(Fu,p,Fd,p,Fnet,p)\n",
    "    axis([0, 700, 0, 1013])\n",
    "    ax = gca()\n",
    "    ax.set_ylim(ax.get_ylim()[::-1])\n",
    "    xlabel('Flux (W m-2)')\n",
    "    ylabel('Pressure (mb)')\n",
    "    title('Longwave Fluxes')\n",
    "    legend(('F_up','F_down','F_net'),loc='best')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OTHER IMPORTANT INFORMATION"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DESCRIPTION OF THE SURFACE REFLECTIVITY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = array([(3250-2600)/2.+2600, (4000-3250)/2.+3250, (4650-4000)/2.+4000, (5150-4600)/2.+4600, (6150-5150)/2.+5150, (7700.-6150)/2.+6150, (8050-7700)/2.+7700, (12850-8050)/2.+8050, (16000-12850)/2.+12850, (22650-16000)/2.+16000, (29000.-22650)/2.+22650, (38000-29000)/2.+29000, (50000-38000)/2.+38000, (2600-800)/2.+800])\n",
    "e = array([0.881, 0.794, 0.738, 0.727, 0.657, 0.660, 0.626, 0.655, 0.791, 0.883, 0.957, 0.958, 0.958, 0.970])\n",
    "figure()\n",
    "plot(10000./v,1-e,'o')\n",
    "xlabel('Wavelength (um)')\n",
    "ylabel('Reflectivity')\n",
    "title('Reflectivity specified for the TRP atmosphere')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DESCRIPTION OF THE SPECTRAL DEPENDENCE OF THE AEROSOL OPTICAL DEPTH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lmbda   = arange(0.,6.,0.1)\n",
    "lmbda1  = 1.0\n",
    "AERPAR1 = 2.184\n",
    "AERPAR2 = 1.00\n",
    "AERPAR3 = 0.00\n",
    "AOD1    = 0.0013     # Example from layer 1.\n",
    "AOD     = AOD1 * (AERPAR2 + AERPAR3 * (lmbda/lmbda1)) / ((AERPAR2 + AERPAR3 - 1) + (lmbda/lmbda1)**AERPAR1)\n",
    "\n",
    "figure()\n",
    "semilogy(lmbda,AOD)\n",
    "xlabel('Wavelength (um)')\n",
    "ylabel('Aerosol Optical Depth')\n",
    "title('Spectral Dependence of the Aerosol Optical Depth')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<i>© Von P. Walden, Washington State University</i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
